---
title: "Distributive Fairness Appeals: Equality, Reciprocity, or Need"
output: 
  html_document:
    toc: false
    toc_depth: 4
    code_folding: hide
    css: style.css
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(rio)
library(dplyr)
library(PerformanceAnalytics)
library(kableExtra)
library(psych)
library(ggplot2)
library(jmv)
library(magick)
library(report)
library(reshape2)
library(ordinal)
library(RVAideMemoire)
library(emmeans)
library(effectsize)
library(repmod)
library(MANOVA.RM)
library(ggrepel)
library(matrixTests)
library(tidyverse)
library(lessR)
library(rstatix)
library(multcomp)
library(jtools)
library(sjPlot)




set.seed(42)

options(scipen=999, digits=3)

paste_na <- function(..., sep = " ") {
  L <- list(...)
  L <- lapply(
    L,
    function(x) {
      x[is.na(x)] <- ""
      x
    }
  )
  out <- gsub(
    paste0("(^", sep, "|", sep, "$)"), "",
    gsub(
      paste0(sep, sep), sep,
      do.call(paste, c(L, list(sep = sep)))
    )
  )
  is.na(out) <- out == ""
  return(out)
}
```

Materials, data, code accompanying the article "Equality, Reciprocity, or Need? Bolstering Welfare Policy Support for Marginalized Groups with Distributive Fairness". 

<br>

# Study 1 {.tabset}
## Data {.tabset}
### Data load
```{r}
s1 <- import("data/Study1.sav")
```

### Data prep
```{r}
s1 <- dplyr::filter(s1, complete.cases(E1A | E1B | E2A | E2B))

s1$SEX <- factor(s1$SEX, levels = c(1, 2), labels = c("Male","Female"))
s1$AGECAT= factor(s1$AGECAT, levels = c(1, 2,3,4,5,6), labels = c("18-24","25-34","35-44", "45-54","55-64","65+"))
s1$EDU= factor(s1$EDU, levels = c(1, 2,3,4), labels = c("Primary","Secondary (no diploma)","Secondary (complete)", "University"))
s1$SIZE= factor(s1$SIZE, levels = c(1, 2,3,4,5), labels = c("less than 1k","1k-4 999","5k-19 999", "20k - 99 999","100k+"))
s1$REG= factor(s1$REG, levels = c(1, 2,3,4,5,6,7,8), labels = c("Bratislavsky","Trnavsky","Trenciansky", "Nitriansky","Zilinsky","Banskobystricky","Presovsky","Kosicky"))

s1$PINCOME <- as.factor(s1$PINCOME)

income <- dplyr::filter(s1, PINCOME != '8' & PINCOME != '9')
income$PINCOME <- as.ordered(income$PINCOME)
median(income$PINCOME)
s1$income <- car::recode(s1$PINCOME,"'1'='below median';'2'='below median'; '8'='NA';'9'='NA'; else = 'above median'")


s1$condition1 <- s1$E1A
s1$condition2 <- s1$E1B
s1$condition3 <- s1$E2A
s1$condition4 <- s1$E2B

s1$condition1[!is.na(s1$condition1)] <- "E1A"
s1$condition2[!is.na(s1$condition2)] <- "E1B"
s1$condition3[!is.na(s1$condition3)] <- "E2A"
s1$condition4[!is.na(s1$condition4)] <- "E2B"

s1$condition <- paste(s1$condition1,s1$condition2,s1$condition3,s1$condition4)
cols <- c("condition1","condition2","condition3","condition4")
s1$condition <- apply(s1[, cols], 1, function(x) toString(na.omit(x)))
s1$condition= factor(s1$condition,  labels = c('S1:Non-Roma settlement', 'S1:Roma settlement','S2:Not including Roma','S2:Including Roma'))

s1$outcome <- paste_na(s1$E1A,s1$E1B,s1$E2A,s1$E2B)
s1$outcome <- as.numeric(s1$outcome)

s1_1 <- dplyr::filter(s1, E1A > 0 | E1B > 0)
s1_1 <- s1_1[, -c(13:14)]
s1_1$group <- s1_1$E1A
s1_1$group[is.na(s1_1$group)] <- "B"
s1_1$group <- car::recode(s1_1$group,"'B'='Roma settlement'; else='Non-Roma settlement'")
s1_1$group <- as.factor(s1_1$group)
s1_1$well <- with(s1_1, pmax(E1A, E1B, na.rm = TRUE))
s1_1$group <- relevel(s1_1$group, "Roma settlement")

s1_2 <- dplyr::filter(s1, E2A > 0 | E2B > 0)
s1_2 <- s1_2[, -c(11:12)]
s1_2$group <- s1_2$E2A
s1_2$group[is.na(s1_2$group)] <- "B"
s1_2$group <- car::recode(s1_2$group,"'B'='including Roma'; else='not including Roma'")
s1_2$policies <- with(s1_2, pmax(E2A, E2B, na.rm = TRUE))
s1_2$group <- as.factor(s1_2$group)
s1_2$group <- relevel(s1_2$group, "including Roma")

#na <- dplyr::filter(s1, !complete.cases(E1A | E1B | E2A | E2B))
```

## Descriptives {.tabset}
### Main outcome 
```{r}
s1_1 %>% dplyr::group_by(group) %>% dplyr::summarise(N = length(well), Min=min(well,na.rm= TRUE),Q1 = quantile(well,probs = .25,na.rm = TRUE),Median = median(well, na.rm = TRUE),Q3 = quantile(well,probs = .75,na.rm = TRUE),Max = max(well,na.rm = TRUE),Mean = mean(well, na.rm = TRUE),SD = sd(well, na.rm = TRUE), Skew = skewness(well, na.rm = TRUE), Kurtosis = kurtosis(well, na.rm = TRUE)) -> s11_table

knitr::kable(s11_table, caption = "Study 1 - 1")%>%
  kable_styling(full_width = F)


s1_2 %>% dplyr::group_by(group) %>% dplyr::summarise(N = length(policies), Min=min(policies,na.rm= TRUE),Q1 = quantile(policies,probs = .25,na.rm = TRUE),Median = median(policies, na.rm = TRUE),Q3 = quantile(policies,probs = .75,na.rm = TRUE),Max = max(policies,na.rm = TRUE),Mean = mean(policies, na.rm = TRUE),SD = sd(policies, na.rm = TRUE), Skew = skewness(policies, na.rm = TRUE), Kurtosis = kurtosis(policies, na.rm = TRUE)) -> s12_table

knitr::kable(s12_table, caption = "Study 1 - 2")%>%
  kable_styling(full_width = F)
```

### Sample characteristics
```{r}
s1_gender <- s1 %>% group_by(SEX) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s1_gender$freq <- round(s1_gender$freq, 3)

s1_gender%>%
kable(caption = "Gender")%>%
  kable_styling(full_width = F)

s1_age <- s1 %>% dplyr::summarise(Min=min(AGE,na.rm= TRUE),Q1 = quantile(AGE,probs = .25,na.rm = TRUE),Median = median(AGE, na.rm = TRUE),Q3 = quantile(AGE,probs = .75,na.rm = TRUE),Max = max(AGE,na.rm = TRUE),Mean = round(mean(AGE, na.rm = TRUE),3),SD = round(sd(AGE, na.rm = TRUE),2),n = n(),Missing = sum(is.na(AGE))) 

s1_age %>%
  kable(caption = "Age") %>%
  kable_styling(full_width = F)

s1_agecat  <- s1 %>% group_by(AGECAT) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s1_agecat$freq <- round(s1_agecat$freq, 3)

s1_agecat%>%
kable(caption = "Age categories")%>%
  kable_styling(full_width = F)

s1_educat <- s1 %>% group_by(EDU) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s1_educat$freq <- round(s1_educat$freq, 3)

s1_educat%>%
kable(caption = "Education categories")%>%
  kable_styling(full_width = F)

s1_reg <- s1 %>% group_by(REG) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s1_reg$freq <- round(s1_reg$freq, 3)

s1_reg%>%
kable(caption = "Region")%>%
  kable_styling(full_width = F)

s1_municip <- s1 %>% group_by(SIZE) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s1_municip$freq <- round(s1_municip$freq, 3)

s1_municip%>%
kable(caption = "Municipality size")%>%
  kable_styling(full_width = F)
```

### Randomization checks
```{r}
s1_check <- dplyr::select(s1, SEX:PINCOME, -AGECAT, condition)
s1_check <- mutate_all(s1_check, function(x) as.numeric(x))

check_s1 <- col_oneway_welch(s1_check[,1:8], s1_check$condition)

check_s1 %>%
  kable(caption = "Randomization checks overall") %>%
  kable_styling(full_width = F)


R1 <- dplyr::select(filter(s1_1, group=="Roma settlement"), SEX:PINCOME, -AGECAT)
R1 <- mutate_all(R1, function(x) as.numeric(x))

nonR1 <- dplyr::select(filter(s1_1, group=="Non-Roma settlement"), SEX:PINCOME, -AGECAT)
nonR1 <- mutate_all(nonR1, function(x) as.numeric(x))

check_S1_1 <- col_t_welch(R1, nonR1)

check_S1_1 %>%
  dplyr::select(mean.x, mean.y,pvalue)  %>%
  kable(caption = "Randomization checks - S1-1") %>%
  kable_styling(full_width = F)


R2 <- dplyr::select(filter(s1_2, group=="including Roma"), SEX:PINCOME, -AGECAT)
R2 <- mutate_all(R1, function(x) as.numeric(x))

nonR2 <- dplyr::select(filter(s1_2, group=="not including Roma"), SEX:PINCOME, -AGECAT)
nonR2 <- mutate_all(nonR2, function(x) as.numeric(x))

check_S1_2 <- col_t_welch(R2, nonR2)

check_S1_2 %>%
  dplyr::select(mean.x, mean.y,pvalue)  %>%
  kable(caption = "Randomization checks - S1-2") %>%
  kable_styling(full_width = F)
```


## Plots {.tabset}

### Settings
```{r}
source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")

raincloud_theme = theme(
text = element_text(size = 11),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 11),
axis.text = element_text(size = 11),
axis.text.x = element_text(vjust = 2),
legend.title=element_text(size=15),
legend.text=element_text(size=15),
legend.position = "right",
plot.title = element_text(lineheight=.7, face="bold", size = 15),
panel.border = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'))
```

### Main outcome
```{r}
descriptives = psych::describeBy(x = s1_1$well, group = s1_1$group)
group = c("Roma settlement","Non-Roma settlement")
means = c(descriptives$'Roma settlement'$mean, descriptives$'Non-Roma settlement'$mean)
se = c(descriptives$'Roma settlement'$se, descriptives$'Non-Roma settlement'$se)
plotdat=data.frame(group, means, se)

#fig <- image_graph(width = 1400, height = 650, res = 192)

(s1_1_main <- ggplot(data = s1_1, aes(y = well, x = group, fill = group)) + geom_flat_violin(data = s1_1, position = position_nudge(x = .2, y = 0), alpha = .5) + geom_point(data = s1_1, aes(y = well, color = group), position = position_jitter(width = .15,height = 0), size = .5, alpha = 0.5) + geom_point(data = plotdat, aes(x = group, y = means), position = position_nudge(x = 0.35), size = 1,inherit.aes = FALSE ) + geom_errorbar(inherit.aes=FALSE, data = plotdat, aes(x = group, ymax= means+(1.96*se), ymin=means+(-1.96*se), y=means), position = position_nudge(x = 0.35), width = 0.3) + expand_limits(x = 5) + guides(fill = FALSE) + guides(color = FALSE) + scale_color_grey() + scale_fill_grey() +
    #scale_color_brewer(palette = "Accent") +
    #scale_fill_brewer(palette = "Accent") +
    theme_bw() + raincloud_theme + labs(y = "How much money should go \nto building a well") + scale_x_discrete(expand = c(0.15,0.2))+ scale_y_continuous(breaks=c(1,3,5,7,9,11)))

# dev.off()  
# out <- image_crop(fig, "750x750+0")
# print(out)

descriptives = psych::describeBy(x = s1_2$policies, group = s1_2$group)
group = c("including Roma","not including Roma")
means = c(descriptives$'including Roma'$mean, descriptives$'not including Roma'$mean)
se = c(descriptives$'including Roma'$se, descriptives$'not including Roma'$se)
plotdat=data.frame(group, means, se)

#fig <- image_graph(width = 1400, height = 650, res = 192)

(s1_2_main <- ggplot(data = s1_2, aes(y = policies, x = group, fill = group)) + geom_flat_violin(data = s1_2, position = position_nudge(x = .2, y = 0), alpha = .5) + geom_point(data = s1_2, aes(y = policies, color = group), position = position_jitter(width = .15,height = 0), size = .5, alpha = 0.5) + geom_point(data = plotdat, aes(x = group, y = means), position = position_nudge(x = 0.35), size = 1,inherit.aes = FALSE ) + geom_errorbar(inherit.aes=FALSE, data = plotdat, aes(x = group, ymax= means+(1.96*se), ymin=means+(-1.96*se), y=means), position = position_nudge(x = 0.35), width = 0.3) + expand_limits(x = 5) + guides(fill = FALSE) + guides(color = FALSE) + scale_color_grey() + scale_fill_grey() + theme_bw() + raincloud_theme + labs(y = "How much money should go\n to policies") + scale_x_discrete(expand = c(0.15,0.2))+  scale_x_discrete(labels = c('Including Roma\n settlements','Not including Roma\n settlements')) + scale_y_continuous(breaks=c(1,3,5,7,9,11)) )

# dev.off()  
# out <- image_crop(fig, "750x750+0")
# print(out)
```

## Main analysis {.tabset}

<b>Study 1 - opinion about the well</b>

```{r, echo=FALSE,results="asis"}
ttest11 <- jmv::ttestIS(data=s1_1, 'well', 'group', welchs = T, meanDiff=T, effectSize=T, mann = T)
report::report(t.test(s1_1$well ~ s1_1$group, paired = F))

cat("<pre>")
print(ttest11)
cat("</pre>")

s1lm_well <- lm(well ~ group, data = s1_1)
#report::report(s1lm_well)
report(s1lm_well)

```

```{r}
s11.ancova <- aov(well ~ group + SEX + AGECAT + EDU + SIZE + REG + PINCOME, data=s1_1)
s11.ancova <-car::Anova(s11.ancova, type="III")
s11eta <- eta_squared(s11.ancova)

knitr::kable(s11.ancova, caption = " ")%>%
  kable_styling(full_width = F)

knitr::kable(s11eta, caption = " ")%>%
  kable_styling(full_width = F)

s1lm_well_cov <- lm(well ~ group + SEX + AGECAT + EDU + SIZE + REG + PINCOME, data = s1_1)
#report::report(s1lm_well_cov)
s11covtab <- report(s1lm_well_cov)
s11covtab<- as.data.frame(s11covtab[["output"]])
#write.csv(s11covtab,"s11covtab.csv")

```


<b>Study 1 - opinions about policies</b>

```{r, echo=FALSE,results="asis"}
ttest12 <- jmv::ttestIS(data=s1_2, 'policies', 'group', welchs = T, meanDiff=T, effectSize=T, mann = T)
report::report(t.test(s1_2$policies ~ s1_2$group, paired = F))

cat("<pre>")
print(ttest12)
cat("</pre>")

s1lm_policies <- lm(policies ~ group, data = s1_2)
report::report(s1lm_policies)
report(s1lm_policies)
```

```{r}
s12.ancova <- aov(policies ~ group + SEX + AGECAT + EDU + SIZE + REG + PINCOME, data=s1_2)
s12.ancova <-car::Anova(s12.ancova, type="III")
s12eta <- eta_squared(s12.ancova)

knitr::kable(s12.ancova, caption = " ")%>%
  kable_styling(full_width = F)

knitr::kable(s12eta, caption = " ")%>%
  kable_styling(full_width = F)

s1lm_policies_cov <- lm(policies ~ group+ SEX + AGECAT + EDU + SIZE + REG + PINCOME, data = s1_2)
#report::report(s1lm_policies_cov)
s12covtab <- report(s1lm_policies_cov)
s12covtab<- as.data.frame(s12covtab[["output"]])
#write.csv(s12covtab,"s12covtab.csv")
```


Interaction with income

```{r}
incomes11 <- dplyr::filter(s1_1, income != 'NA')
lms11 <- lm(well ~ group * income, data = incomes11)
report(lms11)
plot_model(lms11, type = "pred", terms = c("group", "income")) + theme_classic()
#summ <- summ(lms11 )
#summ <- as.data.frame(summ[["coeftable"]])


incomes12 <- dplyr::filter(s1_2, income != 'NA')
lms12 <- lm(policies ~ group * income, data = incomes12)
report(lms12)
plot_model(lms12, type = "pred", terms = c("group", "income")) + theme_classic()
#summ <- summ(lms12 )
#summ <- as.data.frame(summ[["coeftable"]])
```


# Study 2 {.tabset}
## Data {.tabset}
### Data load
```{r}

# s2 <- s2 %>%
#    dplyr::select(-StartDate,-EndDate,-Status,-IPAddress,-Progress,-Duration__in_seconds_,-Finished,-RecordedDate,-ResponseId,-RecipientLastName,-RecipientFirstName,-RecipientEmail,-ExternalReference,-DistributionChannel,-UserLanguage)
# 
# export(s2, "Study2.sav")

s2 <- import("data/Study2.sav")


```

### Data prep
```{r, message=FALSE, results='hide'}
s2$home= factor(s2$home, levels = c(1, 2,3,4), labels = c("home-answered","away-notanswered","home-return", "away-return"))
s2$consent= factor(s2$consent, levels = c(1, 2), labels = c("yes","no"))
s2$etnic_ascr= factor(s2$etnic_ascr, levels = c(1, 2), labels = c("slovak","roma"))
s2$etnicita= factor(s2$etnicita, levels = c(1, 2,3), labels = c("slovak","roma","other"))
s2$age <- as.numeric(s2$age)
s2 <- lessR::recode(age, new_vars="roky_", old=9:72, new=1996:1933, data=s2)
s2$years <- 2019 - s2$roky
s2$gender <- factor(s2$gender, levels = c(1, 2), labels = c("Male","Female"))
s2_consent <- dplyr::filter(s2, consent=="yes")
s2_consent <- dplyr::filter(s2_consent, etnic_ascr != "NA")
s2_consent <- dplyr::filter(s2_consent, eurofondy_iv<5 & skolka_iv<5 & skolka_agree<5 & skolka_vote<5 & skolka_norms < 5 & control < 5 & suma < 5 & praca < 5 & potreba<5 & pila <5 )
```

## Descriptives {.tabset}
### Main outcomes {.tabset}
```{r}
s2_consent %>% dplyr::group_by(etnic_ascr) %>% dplyr::summarise(N = length(control), Min=min(control,na.rm= TRUE),Q1 = quantile(control,probs = .25,na.rm = TRUE),Median = median(control, na.rm = TRUE),Q3 = quantile(control,probs = .75,na.rm = TRUE),Max = max(control,na.rm = TRUE),Mean = mean(control, na.rm = TRUE),SD = sd(control, na.rm = TRUE), Skew = skewness(control, na.rm = TRUE), Kurtosis = kurtosis(control, na.rm = TRUE)) -> s2_control

knitr::kable(s2_control, caption = "Study 2 - control")%>%
  kable_styling(full_width = F)


s2_consent %>% dplyr::group_by(etnic_ascr) %>% dplyr::summarise(N = length(suma), Min=min(suma,na.rm= TRUE),Q1 = quantile(suma,probs = .25,na.rm = TRUE),Median = median(suma, na.rm = TRUE),Q3 = quantile(suma,probs = .75,na.rm = TRUE),Max = max(suma,na.rm = TRUE),Mean = mean(suma, na.rm = TRUE),SD = sd(suma, na.rm = TRUE), Skew = skewness(suma, na.rm = TRUE), Kurtosis = kurtosis(suma, na.rm = TRUE)) -> s2_equal

knitr::kable(s2_equal, caption = "Study 2 - equality")%>%
  kable_styling(full_width = F)

s2_consent %>% dplyr::group_by(etnic_ascr) %>% dplyr::summarise(N = length(praca), Min=min(praca,na.rm= TRUE),Q1 = quantile(praca,probs = .25,na.rm = TRUE),Median = median(praca, na.rm = TRUE),Q3 = quantile(praca,probs = .75,na.rm = TRUE),Max = max(praca,na.rm = TRUE),Mean = mean(praca, na.rm = TRUE),SD = sd(praca, na.rm = TRUE), Skew = skewness(praca, na.rm = TRUE), Kurtosis = kurtosis(praca, na.rm = TRUE)) -> s2_prop

knitr::kable(s2_prop, caption = "Study 2 - proportionality")%>%
  kable_styling(full_width = F)

s2_consent %>% dplyr::group_by(etnic_ascr) %>% dplyr::summarise(N = length(potreba), Min=min(potreba,na.rm= TRUE),Q1 = quantile(potreba,probs = .25,na.rm = TRUE),Median = median(potreba, na.rm = TRUE),Q3 = quantile(potreba,probs = .75,na.rm = TRUE),Max = max(potreba,na.rm = TRUE),Mean = mean(potreba, na.rm = TRUE),SD = sd(potreba, na.rm = TRUE), Skew = skewness(potreba, na.rm = TRUE), Kurtosis = kurtosis(potreba, na.rm = TRUE)) -> s2_need

knitr::kable(s2_need, caption = "Study 2 - need")%>%
  kable_styling(full_width = F)

```

```{r}

s2_consent <- mutate(s2_consent, ID = row_number())

s2_mains <- melt(s2_consent,
        # ID variables - all the variables to keep but not split apart on
    id.vars=c("etnic_ascr", "ID","eurofondy_iv","skolka_iv" ),
        # The source columns
    measure.vars=c("control", "suma", "praca","potreba","pila"),
        # Name of the destination column that will identify the original
        # column that the measurement came from
    variable.name="condition",
    value.name="measurement"
)

s2_mains <- filter(s2_mains, condition!="pila") 

s2_dvs <- s2_mains %>%
  group_by(etnic_ascr,condition)%>%
  dplyr::select(-ID:-skolka_iv) %>%
  get_summary_stats(type = c("mean_sd"))%>%
  dplyr::select(-variable)%>%
  dplyr::mutate(condition = dplyr::recode(condition,
    "control" = "Control",
    "suma" = "Equality",
    "praca" = "Reciprocity",
    "potreba" = "Need"))%>%
  dplyr::mutate(sd = round(sd,2))
      

dvstt <- s2_mains %>% 
#  pivot_longer(contains("att_"))%>% 
    group_by(etnic_ascr)%>% 
    t_test(measurement~condition, detailed = T)      

```


### Sample characteristics
```{r}
s2_consent_gender <- s2_consent %>% group_by(gender) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s2_consent_gender$freq <- round(s2_consent_gender$freq, 3)

s2_consent_gender%>%
kable(caption = "Gender")%>%
  kable_styling(full_width = F)

s2_consent_age <- s2_consent %>% dplyr::summarise(Min=min(years,na.rm= TRUE),Q1 = quantile(years,probs = .25,na.rm = TRUE),Median = median(years, na.rm = TRUE),Q3 = quantile(years,probs = .75,na.rm = TRUE),Max = max(years,na.rm = TRUE),Mean = round(mean(years, na.rm = TRUE),3),SD = round(sd(years, na.rm = TRUE),2),n = n(),Missing = sum(is.na(years))) 

s2_consent_age %>%
  kable(caption = "Age") %>%
  kable_styling(full_width = F)

s2_consent_etnic <- s2_consent %>% group_by(etnic_ascr) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s2_consent_etnic$freq <- round(s2_consent_etnic$freq, 3)

s2_consent_etnic%>%
kable(caption = "Ascribed ethnicity")%>%
  kable_styling(full_width = F)

s2_consent_etnic_self <- s2_consent %>% group_by(etnicita) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s2_consent_etnic_self$freq <- round(s2_consent_etnic_self$freq, 3)

s2_consent_etnic_self%>%
kable(caption = "Self-reported ethnicity")%>%
  kable_styling(full_width = F)
```

## Plots {.tabset}
### Main outcomes
```{r}


#fig <- image_graph(width = 1400, height = 650, res = 192)
    
ggplot(data = s2_mains, aes(y = measurement, x = condition, fill = etnic_ascr)) + geom_flat_violin(data = s2_mains, position = position_nudge(x = .2, y = 0), alpha = .8) + geom_point(data = s2_mains, aes(y = measurement, color = etnic_ascr), position = position_jitter(width = .15,height = .15), size = .5, alpha = 1.2) +
    scale_color_grey() +
    scale_fill_grey(labels=c("Slovak","Roma")) +
    theme_bw() + raincloud_theme + labs(y = "Agreement with building social housing") + scale_x_discrete(expand = c(0.15,0.2),labels=c("Control","Equality","Reciprocity","Need"))+ labs(fill='Ascribed \nethnicity') + guides(color = FALSE) 

#dev.off()  
#out <- image_crop(fig, "1380x750+0")
#print(out)
```


## Main analysis {.tabset}
```{r}
s2_mains$condition <- dplyr::recode(s2_mains$condition, suma = "equal", praca = "reciprocity", potreba="need", control = "control")
s2_mains$measurement <- as.factor(s2_mains$measurement)

model <- clmm(measurement ~ condition * etnic_ascr + (1|ID),data = s2_mains, threshold = "equidistant", Hess = T)
summary(model)
Anova.clmm(model, type = "II")

emmip <- emmip(model, condition ~ etnic_ascr, CIs = T, engine = "ggplot",style = "factor", plotit = F, mode= "mean.class")

names(emmip)[8] <- "Condition"
emmip$Condition <- dplyr::recode(emmip$Condition , 'control' = "Control", "reciprocity" = 'Reciprocity', "need" = 'Need', "equal" = 'Equality')
emmip$xvar <- dplyr::recode(emmip$xvar , 'slovak' = "Slovak", "roma" = 'Roma')

(emmplot <- ggplot(data = emmip, aes(y= yvar, x = xvar, group = Condition)) + geom_point(size=1.3, position=position_dodge(width=0.1), aes(color = Condition)) + #geom_line(size=1.3,position=position_dodge(width=0.1),aes(linetype=Condition, color=Condition))+  
    geom_linerange(size=1.3,aes(x = xvar,
                     ymin = LCL,
                     ymax = UCL, color = Condition),
                 show.legend = FALSE,position=position_dodge(width=0.1))+
    scale_color_grey(start = 0.2, end = 0.7) +
    scale_fill_grey(start = 0.2, end = 0.7)+theme_classic()  + labs(y = "Predicted average agreement", x = 'Ascribed ethnicity') +
  geom_text_repel(aes(label = Condition, colour = Condition),
    data          = subset(emmip, xvar == "Slovak"),
    nudge_x       = -0.25,
    segment.size  = 0.2,
    direction     = "y") +
  geom_text_repel(aes(label = Condition, colour = Condition),
    data          = subset(emmip, xvar == "Roma"),
    nudge_x       = 0.25,
    segment.size  = 0.2,
    direction     = "y") + theme(legend.position="none")+
  scale_linetype_manual(values=c("twodash", "dotted","longdash", "solid" ))+
    theme(text = element_text(size=15)) + ylim(1,4) )

#ggsave("plot.png", width = 6, height = 5)

emmip%>%
kable(caption = "Predicted average agreement")%>%
  kable_styling(full_width = F)
```


# Study 3 {.tabset}
## Data {.tabset}
### Data load
```{r}
s3 <- import("data/Study3.sav")
```

### Data prep
```{r, message=FALSE, results='hide'}
s3$group= factor(s3$SKUPINA, levels = c(1, 2,3,4), labels = c("Control","Equality", "Proportionality", "Need"))
s3$SEX <- factor(s3$SEX, levels = c(1, 2), labels = c("Male","Female"))
s3$AGECAT= factor(s3$AGECAT, levels = c(1,2,3,4,5,6), labels = c("18-24","25-34","35-44", "45-54","55-64","65+"))
s3$EDU= factor(s3$EDU, levels = c(1, 2,3,4), labels = c("Primary","Secondary (no diploma)","Secondary (complete)", "University"))
s3$SIZE= factor(s3$SIZE, levels = c(1, 2,3,4,5), labels = c("less than 1k","1k-4 999","5k-19 999", "20k - 99 999","100k+"))
s3$REG= factor(s3$REG, levels = c(1, 2,3,4,5,6,7,8), labels = c("Bratislavsky","Trnavsky","Trenciansky", "Nitriansky","Zilinsky","Banskobystricky","Presovsky","Kosicky"))


s3$income <- car::recode(s3$PINCOME,"'1'='below median';'2'='below median';'8'='NA';'9'='NA'; else = 'above median'")
```

## Descriptives {.tabset}
### Main outcomes {.tabset}
```{r}
s3 %>% dplyr::group_by(group) %>% dplyr::summarise(N = length(E3), Min=min(E3,na.rm= TRUE),Q1 = quantile(E3,probs = .25,na.rm = TRUE),Median = median(E3, na.rm = TRUE),Q3 = quantile(E3,probs = .75,na.rm = TRUE),Max = max(E3,na.rm = TRUE),Mean = mean(E3, na.rm = TRUE),SD = sd(E3, na.rm = TRUE), Skew = skewness(E3, na.rm = TRUE), Kurtosis = kurtosis(E3, na.rm = TRUE)) -> s3_personal

knitr::kable(s3_personal, caption = "Study 3 - personal agreement")%>%
  kable_styling(full_width = F)

s3 %>% dplyr::group_by(group) %>% dplyr::summarise(N = length(E4), Min=min(E4,na.rm= TRUE),Q1 = quantile(E4,probs = .25,na.rm = TRUE),Median = median(E4, na.rm = TRUE),Q3 = quantile(E4,probs = .75,na.rm = TRUE),Max = max(E4,na.rm = TRUE),Mean = mean(E4, na.rm = TRUE),SD = sd(E4, na.rm = TRUE), Skew = skewness(E4, na.rm = TRUE), Kurtosis = kurtosis(E4, na.rm = TRUE)) -> s3_majority

knitr::kable(s3_majority, caption = "Study 3 - agreement of majority")%>%
  kable_styling(full_width = F)
```

### Sample characteristics
```{r}
s3_gender <- s3 %>% group_by(SEX) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s3_gender$freq <- round(s3_gender$freq, 3)

s3_gender%>%
kable(caption = "Gender")%>%
  kable_styling(full_width = F)

s3_age <- s3 %>% dplyr::summarise(Min=min(AGE,na.rm= TRUE),Q1 = quantile(AGE,probs = .25,na.rm = TRUE),Median = median(AGE, na.rm = TRUE),Q3 = quantile(AGE,probs = .75,na.rm = TRUE),Max = max(AGE,na.rm = TRUE),Mean = round(mean(AGE, na.rm = TRUE),3),SD = round(sd(AGE, na.rm = TRUE),2),n = n(),Missing = sum(is.na(AGE))) 

s3_age %>%
  kable(caption = "Age") %>%
  kable_styling(full_width = F)

s3_agecat  <- s3 %>% group_by(AGECAT) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s3_agecat$freq <- round(s3_agecat$freq, 3)

s3_agecat%>%
kable(caption = "Age categories")%>%
  kable_styling(full_width = F)

s3_educat <- s3 %>% group_by(EDU) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s3_educat$freq <- round(s3_educat$freq, 3)

s3_educat%>%
kable(caption = "Education categories")%>%
  kable_styling(full_width = F)

s3_reg <- s3 %>% group_by(REG) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s3_reg$freq <- round(s3_reg$freq, 3)

s3_reg%>%
kable(caption = "Region")%>%
  kable_styling(full_width = F)

s3_municip <- s3 %>% group_by(SIZE) %>% 
dplyr::summarise(n = n()) %>%
  mutate(freq = n / sum(n))
s3_municip$freq <- round(s3_municip$freq, 3)

s3_municip%>%
kable(caption = "Municipality size")%>%
  kable_styling(full_width = F)
```

### Randomization checks
```{r}
s3_check <- dplyr::select(s3, SEX:PINCOME, -AGECAT, group)
s3_check <- mutate_all(s3_check, function(x) as.numeric(x))

s3_check <- col_oneway_welch(s3_check[,1:5], s3_check$group)

s3_check %>%
  kable(caption = "Randomization checks overall") %>%
  kable_styling(full_width = F)

table(s3$SEX, s3$group)
```


## Plots {.tabset}
### Main outcomes
```{r}
descriptives = psych::describeBy(x = s3$E3, group = s3$group)
group = c("Control","Equality", "Proportionality", "Need")
means = c(descriptives$'Control'$mean, descriptives$'Equality'$mean,descriptives$'Proportionality'$mean, descriptives$'Need'$mean)
se = c(descriptives$'Control'$se, descriptives$'Equality'$se,descriptives$'Proportionality'$se, descriptives$'Need'$se)
plotdat=data.frame(group, means, se)

(s3_personal <- ggplot(data = s3, aes(y = E3, x = group, fill = group)) + geom_flat_violin(data = s3, position = position_nudge(x = .2, y = 0), alpha = .5) + geom_point(data = s3, aes(y = E3, color = group), position = position_jitter(width = .15,height = .15), size = .5, alpha = 0.5) + geom_point(data = plotdat, aes(x = group, y = means), position = position_nudge(x = 0.35), size = 1,inherit.aes = FALSE ) + geom_errorbar(inherit.aes=FALSE, data = plotdat, aes(x = group, ymax= means+(1.96*se), ymin=means+(-1.96*se), y=means), position = position_nudge(x = 0.35), width = 0.3) + expand_limits(x = 5) + guides(fill = FALSE) + guides(color = FALSE) + scale_color_grey() + scale_fill_grey() +
    #scale_color_brewer(palette = "Accent") +
    #scale_fill_brewer(palette = "Accent") +
    theme_bw() + raincloud_theme + labs(y = "Personal agreement with building social housing \n(1 = completely disagree)") + scale_x_discrete(expand = c(0.15,0.2),labels=c("Control","Equality", 'Proportionality'="Reciprocity", "Need")))

#ggsave("plot.png", width = 6, height = 4)

descriptives = psych::describeBy(x = s3$E4, group = s3$group)
group = c("Control","Equality", "Proportionality", "Need")
means = c(descriptives$'Control'$mean, descriptives$'Equality'$mean,descriptives$'Proportionality'$mean, descriptives$'Need'$mean)
se = c(descriptives$'Control'$se, descriptives$'Equality'$se,descriptives$'Proportionality'$se, descriptives$'Need'$se)
plotdat=data.frame(group, means, se)

(s3_majority <- ggplot(data = s3, aes(y = E4, x = group, fill = group)) + geom_flat_violin(data = s3, position = position_nudge(x = .2, y = 0), alpha = .5) + geom_point(data = s3, aes(y = E4, color = group), position = position_jitter(width = .15,height = .15), size = .5, alpha = 0.5) + geom_point(data = plotdat, aes(x = group, y = means), position = position_nudge(x = 0.35), size = 1,inherit.aes = FALSE ) + geom_errorbar(inherit.aes=FALSE, data = plotdat, aes(x = group, ymax= means+(1.96*se), ymin=means+(-1.96*se), y=means), position = position_nudge(x = 0.35), width = 0.3) + expand_limits(x = 5) + guides(fill = FALSE) + guides(color = FALSE) + scale_color_grey() + scale_fill_grey() +
    #scale_color_brewer(palette = "Accent") +
    #scale_fill_brewer(palette = "Accent") +
    theme_bw() + raincloud_theme + labs(y = " Majority agreement with building social housing \n(1 = completely disagree)") + scale_x_discrete(expand = c(0.15,0.2),labels=c("Control","Equality", 'Proportionality'="Reciprocity", "Need")))

#ggsave("plot.png", width = 6, height = 4)
```

## Main analysis - separate ordinal regression {.tabset}

```{r}

s3$E1 <- as.ordered(s3$E1)
s3$E2 <- as.ordered(s3$E2)
s3$E3 <- as.ordered(s3$E3)
s3$E4 <- as.ordered(s3$E4)


s3personal  <- polr(E3 ~ group, data = s3, Hess=TRUE)
summary(s3personal)
ctable <- coef(summary(s3personal))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = p))
ci <- confint(s3personal)
exp(cbind(OR = coef(s3personal), ci))

s3majority  <- polr(E4 ~ group, data = s3, Hess=TRUE)
summary(s3majority)
ctable <- coef(summary(s3majority))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = p))
ci <- confint(s3majority)
exp(cbind(OR = coef(s3majority), ci))


# Re-level to have Proportionality as baseline

s3$group <- relevel(s3$group, ref = "Proportionality")
s3personal  <- polr(E3 ~ group, data = s3, Hess=TRUE)
summary(s3personal)
ctable <- coef(summary(s3personal))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = p))
ci <- confint(s3personal)
exp(cbind(OR = coef(s3personal), ci))


s3majority  <- polr(E4 ~ group, data = s3, Hess=TRUE)
summary(s3majority)
ctable <- coef(summary(s3majority))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = p))
ci <- confint(s3majority)
exp(cbind(OR = coef(s3majority), ci))

```


## Main analysis - ordinal regression, interactions with income {.tabset}

```{r}
s3income <- dplyr::filter(s3, income != 'NA')

s3income$group <- relevel(s3income$group, ref = "Control")

s3income3  <- polr(E3 ~ group * income, data = s3income, Hess=TRUE)
summary(s3income3)
ctable <- coef(summary(s3income3))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = p))
ci <- confint(s3income3)
exp(cbind(OR = coef(s3income3), ci))


s3income4  <- polr(E4 ~ group * income, data = s3income, Hess=TRUE)
summary(s3income4)
ctable <- coef(summary(s3income4))
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
(ctable <- cbind(ctable, "p value" = p))
ci <- confint(s3income4)
exp(cbind(OR = coef(s3income4), ci))

```



Interaction with EU fund attitudes and perceived competitiveness
```{r}
s3$E1_coll <- car::recode(s3$E1,"'1'='2'")
s3$E2_coll <- car::recode(s3$E2,"'1'='2'")

s3$E1_coll <- as.ordered(s3$E1_coll)
s3$E2_coll <- as.ordered(s3$E2_coll)

s3_fit_all <- MANOVA.wide(cbind(E3, E4) ~ group * E1_coll * E2_coll, data = s3, iter = 1000, CPU = 1)

summary(s3_fit_all)

simCI(s3_fit_all, contrast = "pairwise", type = "Tukey", interaction = F, factor ="group")

```


# Session Info
```{r}
sessionInfo()
```
